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Abstract. In this paper we present and test a full discretization of all elements of the Calderon 

Calculus (layer potentials and integral operators) for the Helmholtz equation in smooth closed curves 

in the plane. The resulting integral equations provide approximations of order three for all variables 

involved. Test are shown for a wide array of direct, indirect and combined field integral equation at 

^f^ fixed frequency and for a Convolution Quadrature based approximation in the time domain. 
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1. Introduction. This paper introduces and tests a fully discrete Calderon Cal- 
culus for two dimensional acoustic waves, time-harmonic and transient. We present a 
^s^ novel simultaneous discretization of the two layer potentials and four integral opera- 

tors associated to the Helmholtz equation at fixed frequency on a collection of smooth 
' ^ ' non-intersecting parametrizable closed curves in the plane. The method's vocation is 

^N simplicity, and we can assert with some confidence, that it will not be easy to find 

^_j another instance of such a simple method of reasonable order (order three in all vari- 

fH ablcs, in strong norms), with so little computational and programming requirements. 

We do not make any claims, though, on the ability of this method to work on prob- 
lems at high frequencies, and we are by no means competitors of sophisticated high 
order methods that look for fine details in complicated geometries. We do, however, 
claim that the set of tools exposed in this paper works for many other integral op- 
erators (experiments on the Laplace equation have been carried out by the authors 
^ as a prototyping tool), and we are working on the extension of these ideas to some 

more general problems. It is important to emphasize that we are not discretizing a 
— 1 particular integral equation, so we are not worrying on whether one formulation is 

r^^ better than another, or whether there are resonances. As we show in the examples, 

./ the discrete operators and potentials can be used to build any of the best known 

f~^ direct, indirect, and combined field integral equations for exterior problems, as well 

fT^ as more complicated systems of integral equations for transmission problems. The 

T-H time domain extension is carried out by using a variable complex frequency and the 

►^ Convolution Quadrature technology of Christian Lubich |5H |5| . 

The method works in a relatively simple way. On a parametric curve, several sets 
of points and normal vectors are sampled. They are harvested using three staggered 
^ uniform grids in parameter space. One grid is used to create sources and two grids 

are used for simultaneous (averaged) observation. Once each curve is sampled at 
the discrete level, merging information to create a unified discrete set is an easy 
task. The second step is the automatic creation of potentials and operators using 
direct evaluations of the kernel functions: all numerical integration processes are done 
explicitly in the method, and all equations and right-hand sides are fully discrete. 
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The ideas behind these methods go back to a very simple quadrature method of 
order two by Saranen and Schroderus US] , later generalized [5] and improved to third 
order of convergence |T3j . The treatment of the associated hypersingular operator is 
surprisingly simple as well: using the integration by parts formula that is common to 
Galerkin discretizations of the hypersingular integral equation, the paper [12] found a 
collection of fully discrete discretizations of order one and two for this hypersingular 
equation. Some additional work allowed us to put together the first Calderon Calculus 
of order two in [TT] . This set of discrete operators is heavily asymmetric and has the 
disadvantage of requiring sampling of second derivatives of the parametrization of the 
curve due to the evaluation of the double layer operator on its diagonal. The current 
set of methods mixes the discoveries of [TT] and [13] to create a discrete set of order 
three that is even simpler than the order two collection. 

The discrete set is, as a matter of fact, a Nystrom (quadrature) discretization 
of the integral operators, avoiding evaluation of singular kernels on the diagonal, 
looking for superconvergent location of observation points, and mixing observation 
grids to partially symmetrize the method, and achieve order three. However, the 
method can be better understood as a full discretization, with carefully chosen low 
order numerical integration, of a non- conforming Petrov-Galerkin discretization of 
the integral operators. The methods will be tested on a wide set of integral equations 
for exterior and transmission problems, and on a time-domain scattering problem. 
We will also test condition numbers of the different formulations and the possibility 
of using Calderon preconditioning. 

Some discussion on the literature. Nystrom methods [M1[T] are the most popular 
choices for integral equations of the second kind. For integral equations of the second 
kind with smooth periodic kernels, the trapezoidal rule gives rise to a very powerful 
method which converges superalgebraically [TU Chapter 12]. Periodic weakly singular 
integral equations of the second kind (as those that arise from the Helmholtz equation 
on smooth parametrizable domains in the plane) are also amenable to simple methods 
with superalgebraic or exponential order of convergence [5] Section 3.5] (see also 
[201 US]). A comparison of Nystrom methods in the plane has been carried out in [IB] . 
The three dimensional case is much more involved and consequently less developed. 
There are Nystrom schemes for equations with weakly singular kernels, like those 
of Bruno and Kunyaski [H |3] and Wienert [H EH [H]. Finally, the QBX methods 
[TBI dU [TB] , originally designed to compute layer potentials close to the boundary, 
are proving to be useful tools to create Nystrom methods for weakly singular integral 
equations in two and three dimensions. 

2. Parametrized Calderon Calculus. Let x : M ^ F C M^ be a smooth 
(x e C°°(M)) regular (|x'(i)| 7^ for all t) 1-periodic (x(t + 1) = x(t) for all t) 
positively oriented parametrization of a simple (x(i) 7^ x(t) if < t < t < 1) closed 
curve in the plane. We consider the parametrized non-normalized normal vector field 
n{t) :— (x2(i), —x'i{t)). The curve F divides the plane into a bounded interior domain 
ri_ and its unbounded exterior $7+ . Given a function [/ : M^ \ F — > C that is smooth 
enough on both sides of the interface F, we will write U^ for its restrictions of r2±. 
Restrictions (traces) and normal derivatives on the boundary will be defined as: 

7±[/:=[/±|rox, 9±C/:-((VC/±)|rox)-n. (2.1) 

As defined, these restrictions to the boundary define periodic functions and that 
the normal derivative is actually a directional derivative with respect to the non- 
normalized outward pointing normal vector field n. 
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Given complex valued sufficiently smooth 1-periodic functions, we can define the 
single and double potentials on T with the formulas 



I '•1 



4./0 



Sv){^):^- Hi,'\k\z^^mv{t)dt, (2.2a) 



(D V')(z) := ^ ^' H['\k\z x(t)|) ^"^/_^'yV w dt. (2.2b) 

These functions are defined for all z S K''\r. It is well known (it actually follows from 
a very simple computation) that for any r], ^ the function U := 877 + 0-0 € C°°(M^ \r) 
solves 

AU + k^U = OmR^\r, lim |z|i/2(vC/(z) • (yir z) - ifc[/(z)) = 0, (2.3) 

|z|— >oo V ' ' / 

that is, [/ is a radiating solution of the Hclmholtz equation. Smoothness of U as we 
approach the interface F depends on smoothness of the densities rj and ip. Reciprocally, 
given a solution of (2.3), we can write 

C/ = Sia„(7]-Dl7C/l, (2.4) 

where the jump operators are defined as 



Uniqueness of the representation (2.4) for the solutions of (2.3) implies the jump 
relations of potentials 

|7S7?l-0, ldnSv]=v, [7D^1=-^, [9„D^1 = 0. (2.5) 

These jump properties motivate the introduction of the four operators on the bound- 
ary T: 

yV-=l7Sv} = 7^STi, Jv-={{dnSv}, (2.6a) 

Ki^:={{7D0}, WV':--{{a„D7A:& = -a±DV, (2.6b) 

where 

= Uj-u + j+u), {{ant/S:-i(a„c/ + 9+[/). 



The operators V and K are the single and double layer operators respectively, J is the 
adjoint double layer operator, and W is the hypersingular operator for the Helmholtz 
equation. The first three of these operators admit integral expressions: 



* 1^ 'r(l) 



(Vr/)(T):=-y^ H^'>ik\^iT)-^it)\Mt)dt, (2.7a) 

imir) := ^ f H[^\k\.ir) - ^WD^^'^^^f^^^m ^t, (2.7b) 
4 Jo |x(T)-x(i)| 

(J.)(r) := ^ /'i/f)(fc|x(.)^x(.)|) (-(;)--(-»."(-^ (t)dt. (2.7c) 
4 Jo |x(t) -x(i)| 



It is clear from here that J = K*. We wiU keep a different notation though, for reasons 
that will become apparent when we discretize them in a non-symmetric form. The 
operator W admits an expression in the form of an integro-differential operator: 



where 



WV' 



-(W/) 



k^Y, 



^, 



(2.7d) 



(V„V^)(r) 



i/«(fc|x(r)-x(t)|)(n(i).n(T))^(i)dt. 



(2.7e) 



The representation formula (2.4) for all radiating solutions of the Helmholtz equation 



(2.3), together with the jump properties of the potentials (2.5) and the definitions of 



the boundary integral operators by averaging (2.6), determines a set of rules (a calcu- 



lus) that generates a diverse collection of representation formulas, potential ansatzs, 
and integral equations associated to the solution of interior, exterior and transmission 
problems for the Helmholtz equation. The following matrices of operators 



7 
5± 



[D -S]=± 



I 
I 



K -V 
-W -J 



(2.8) 



collect the exterior/interior Cauchy values of the layer potentials. They constitute 
the exterior/interior Calderon projectors associated to the Helmholtz equation. The 
systematic use of these potentials and operators to build integral equations will be 
explored in Section |6] At this point, let us emphasize the fact that we are striving for 
a full discretization of the entire set of potentials (2.2) and operators (2.7), including 



also discretization of the restriction operators (2.1 ) that will be needed to sample, at 
the discrete level, incoming incident waves. 

The case of multiple scatterers. Assume that ri,...,rM are pairwise disjoint 
curves, parametrized as above by smooth 1-periodic functions x^. All potentials and 
operators can be easily defined for vectors of densities (771, . . . , 77m) and (i/)!, . . . , iPm)- 
The integral operators then become matrices of integral operators. For instance, we 
have operators of the form 



r(i) 



H^„'>{k\Mr)-^mmVm{t)dt, 



£,m^l,...,M. 



3. Fully discrete method. 



3.1. Geometry. For one single curve parametrized with x as in Section [2] we 
proceed as follows. We take a positive integer N and define h := 1/N. Next we define 
the discrete parameter points tj :— hj and the values 



nij := x(tj), hj :— x{tj — h/2), 



hn{tj), jeZAr:={l,...,iV}. 



The notation m.j and hj makes reference to midpoints and breakpoints of a boundary 
element mesh that is implicit to this method (see Section [4]) . In addition to these 
sampled quantities, we need two index-based functions that provide the next and 
previous index modulo N: the next-index function is n : Z^r — > Zjv given by 



nU) ■= 



j + 1, l<j<iV-l, 
1, J = N, 
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while p := n^^ . Merging geometric information from several curves is easy: after sam- 
pling two curves Fi and r2 with Ni and N2 elements respectively, midpoints, break- 
points, and normals are collected in lists with N = Ni + N2 elements, by appending 
the information of T2 after the information of Fi. We then create the next- index and 
previous-index functions by juxtaposing the two existing functions: 



n{j) 




P 



1 < J < iVi - 1, 

J = Ni, 

Ni + l<j<Ni+N2-l, 

j = Ni+ N2, 

This merging process can be applied to any finite number of curves, each one dis- 
cretized (sampled) with a different number of points. The quantity h appears only 
at the time of collecting information from a particular curve and is incorporated to 
quantities related to first derivatives of the parametrization. However, at the time of 
merging, h is absent from any expression. From this moment on, n : Z^r — >• Zjv is a 
permutation of Z^v and p = n^^ . 

3.2. Discrete potentials. The discrete version of the single and double layer 



potentials (2.2) is defined by using linear combinations of monopoles and dipoles: 



Given two vectors rj — {rji, . 
tials 



and 



ik 



r(i) 



D,{z):^-H['>ik\z-m,\)- 



- mJ 



{^pi,...,ipNy 



^N 



N 



Sft(z)T7:=^77j$j(z), 



the discrete poten- 



(3.1a) 



D/,(z)t/. 



N 

k E V'j(^pO)(^) + 22^.7- (z) + A.o-)(z)) 

N 

h H'^^P(J) + 22^-^ + i/'„(j))£ij(z). 



(3.1b) 



define solutions of (2.3) 



A quadrature related matrix. Let us consider the N x N matrix Q given by 

Qi,t = 15, Qi,n(j) = Qi,p{i) = ^, Qjj = otherwise. (3.2) 

When the geometry proceeds from a single sampled curve, and therefore the next- 
index function is just a right-shift modulo N, Q is the circulant symmetric matrix 

22 1 1 

1 22 1 



Q 



1 

24 



1 



22 
1 



1 
22 



In general, Q is block diagonal with blocks of the above form, one for each of the 
curves. This matrix is related to a quadrature formula that will be introduced in 
Section [4] It is clear that (3.1b) is just a linear combination of dipoles, where either 
the coefficients are premultiplied by the matrix Q, or the dipoles themselves are mixed 
using this matrix. 
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3.3. Observation grids and mixing matrices. Since the integral operators in 
(2.7) have singularities at r = t, we are forced to use a different discrete set for testing. 
We start by defining two sets of discrete samples. For a single curve parametrized 
with X, we use the same N and h = 1/A^ to define 



mi 



x{tj±h/6), hf --xltj -h/2±h/6), nf := hx{tj ± h/6), j eZ 



^N- 



As in Section |3.1[ observations on finite collections of curves are merged in a simple 
way. We demand that the number of discretization and observation points on each 
curve coincides, although it can be taken to be different on different curves. 

Instead of directly averaging values from both possible choices, we will be con- 
sidering a more general mixture of the two grids. We start with the N x N matrix 
P+ = P+(a) with elements 



P+. •= ia 

^ 1,1 ■ 2 ' 



P+,. 



= i(l-a), 



PJ'j =0 otherwise. 



The parameter a > will be discussed in Sectional We also let P^ := (P+)^. For 
the case of a single curve (when n is the right-shift modulo N) , we show two particular 
cases of interest: 



p^(i^ 



1 

12 



P+(l) 



il. 



Given two vectors ^* e 1 

(p+i+ + p-e 



1 5 
]I^, it is easy to see that 

1. = i((i-«)Cw + <^+<' 



Similarly, the i-th element of Q(P+^+ + P"^") = P+Q^^ 



(i-«)C(,))- 

P^Q^~, namely 



(3.3) 



s((i-«)e 



pH'O 



ae 



p{i) 



(22 - 21a)e+,) + (21a + 1)^: 



+{21a + 1)C+ + (22 - 21a)q,) + aC(,) + (1 - «)&(,)) , 



(3.4) 



is a weighted local average of the values around the index i. 

The testing part of the discrete Calderon Calculus is applied upon an incident 
wave. At this point, this is just a function [/'"'^ : M^ ^- C that is smooth around the 
collection of curves, so that we can evaluate 



(m±))^ 



(3.5a) 
(3.5b) 



/3±:=-(C/'-(m±),...,C/'- 

/3± := -(VC/'-(m±) • n±, . . . , VC/'-(m± ) • n± )^ 

We finally define the observation of the incident wave and its normal derivative with 

f3o:=F+f3++p-f3^, A := Q(P+/3+ + p-/3r)- (3.5c) 

3.4. Discrete operators. The discrete operators are defined using the geomet- 
ric elements of Section [3T] in the integration variable and the observation grids of Sec- 
tion [TS] in the test variable. The subscript h will be used to denote discretization. In 
the case of several curves {Fi, . . . , Fa/}, we can consider that h := (l/A^i, . . . , 1/Nm), 
although this is not relevant for the exposition of the methods. 
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Following (2.7) we define two sets of discrete operators (based on the principal 
sampling of the geometry, tested on both ± observation grids). We start with the 
three integral operators 



V 



»j 



:=^ijW(fc|m±-m,|), 



K,^ 



I k 
I k 



H[^\k\raf - m^- 



(m± 



Im,* — m. 



*j 



i/«(fc|m,-m±|) 



(mj 



m, 



m. 



(3.6a) 
(3.6b) 

(3.6c) 



Following (2.7d), the discretization of W separates the discretization of the principal 
part 



W? 



V± 

n{i),n(j) 



V; 



(Oj 



V 



± 

iMj) 



■V 



r(i) 



V?.^=7^r(fc|bf-b,|), (3.6d) 



from the more regular logarithmic term in (2.7e) 



¥=■ 



(nf-n,)V^ 



(3.6e) 



If V^,K^, J^, W^, and V|^^ are the above matrices, we define the matrices of the 
discrete Calderon Calculus by 






-V+ 



p+vr + P"V 



h 

-P+K+Q 
= P+QJ+^ 

= P+W+ 4 
= P+W+ 4 



h ' 

+ P-K-Q = (P+K+ + P-K-)Q, 
-P-QJ^=Q(P+J++P-J^) 
P-W- - fc2(P+QV+,Q + P-QV„_,Q) 

p-w,;-fc2Q(p+v+ +p-v„,)Q. 



(3.7a) 
(3.7b) 
(3.7c) 
(3.7d) 
(3.7e) 



The Calderon projectors (2.8) include the action of two identity operators. Both of 
them will be approximated by the following mass matrix M ~ M(a) 



M, 



Kl + 3a), M,,p(,)=M, 



j^(7-6a), M,,j=0 otherwise. (3.8) 



The simplest method corresponds to a = 1. In this case P* = 5 I ^ind, apart from the 
action of the matrix Q (related to quadrature), we are just averaging sets of equations 
on the two grids. However, even in this simple case, the mass matrix has a circulant 
tridiagonal structure. 

4. From Nystrom to Petrov-Galerkin. In this section we reinterpret all the 
matrices and testing of right-hand sides given in Section [3] as non-conforming Petrov- 
Galerkin method with numerical quadrature. This will be done for the case of a single 
curve, where we are working with a single equation and parametric unit interval (1- 
periodic real line). When there are M curves, M copies of the unit interval have to 
be used. The details just became slightly more cumbersome, but all the following 
arguments can be extended readily. 

Discrete functions and spaces. We start by setting some notation. Given z S M, 
we write Sz to denote the 1-periodic Dirac delta distribution at z, that is, the Dirac 
comb supported on z + Z. Given an open interval /, of length less than one, we write 
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G ^- — O 

O >' O &^ 



ffi N: 



Figure 4.1. The shape of the combination of Dirac deltas S* and the piecewise constant function 
Xi ■ The plot is given for the choice a = 5/6. A piecewise linear function and a quadratic spline are 
shown in the background. They are at the origin of the choice of coefficients for the distribution 5^ . 



Xi to denote the 1-periodic function that coincides with the characteristic function of 
/ on a unit length interval containing /. We then write 



St., sf 



^ti±h/6, Xi '■— X{ti-h/2,ti+h/2), Xi 



X(ti±h/e-h/2,ti±h/6+h/2)- 



Next we define the Dirac fork (see (3.3 1 to recognize the corresponding coefficients) 

(4.1) 



S* := i ((1 - a)S+_, + aS^ + a5+ + (1 - a)5,+i) , 
and the ziggurat-shaped piecewise constant functions 

X*t := 5 ((1 ~ Oi)xi-i + ax7 + ^^xt + (1 " a)Xi+i) • 



(4.2) 



Figure 4.1 shows the shapes of the basic test functions for the particular case a — 5/6. 
Using momentarily the notation s^ := ti — hjl ± /i/6, it is easy to note that, in the 
sense of periodic distributions, 



iixt = s,f-K_ 



± 

i + 1 



and 



Ax* = i((l-a)<5,+_^ + 



a8- + a8^ 



+ (1-«)'5.^:;J 



-if(l-a)(5 + +a(5 - + a5 + + (1 - a)(5 - V 



This shows how, in the same way that characteristic functions arise from integrating 
two consecutive deltas with opposite signs, the ziggurat functions arise from integrat- 
ing Dirac forks. Four spaces are relevant for what follows: 



T/i :== span{(5i : i € In}-, 
Sh ■■= spanjxi : i G Zjy}, 



T,; := span{(5.* : i £ Zn}, (4.3a) 

5*;^ := spanjx- : i e In)- (4.3b) 



Note that Sh is just the space of periodic piecewise constant functions on a uniform 
mesh with mesh-size h and {ti\ as midpoints of the mesh elements. The T spaces will 
be non-conforming discretizations of H~^^^ Sobolev spaces, while the S spaces are 
non-conforming approximations of if ^/^. The • spaces will do the job of test spaces, 
while the unscripted spaces will be the trial spaces. 

Interactions of deltas and characteristic functions. We define the actions of deltas 
with characteristic functions with the formulas 



^,S.) 



{X^+1,S^) 



ix^ 
= {xti,s,) 

{S^,xf) 



J_ 
12 



{S^,X.), 



= %-a=: ((5,ti,X*) = (<5,+i,Xj), 



= 



=:('5f,X.) 



otherwise. 



(4.4a) 
(4.4b) 
(4.4c) 



The otherwise case above has to be understood modulo N. This interaction will be 
explained in Section [SJ It is clear that the first line of ( 4.4 ) enforces the second, if 
we want some kind of consistency of our formulas with respect to translations in the 



origin of the real line. The interactions (4.4| and the definitions (4.1), (4.2| imply 
that (see ^M) 



{5lXj) ^ {Xl5,) 



= |(l + 3a) = M,,„ 
= 0, otherwise. 



In other words, the matrix M is the matrix that represents the 'dualities' 5^ x Th and 
T^ X Sh if (4.4) is imposed. It is to be noticed that in the simplest case (a = 1), the 
interaction of a Dirac delta with a characteristic function is forced to be negative on 
neighboring elements ( 4.4b ) . We will discuss these choices in Section [s] 

First collection of discrete elements. While the angled bracket (linear in both 
components) is used for the concrete interactions of piecewise constant functions and 
Dirac delta distributions, from now on we will use curly brackets (linear in both 
components as well) for the following situations 



{5,,4>}:^^{z), {xi,4>}-= j m<^t- 



This can be applied as long as the right-hand side of the expression is meaningful. 
We can then define the following bilinear forms 



K^Tu3{iil,ilh) 


^-> {m^Vt^/i}, 


(4.5a) 


SlxSH3{(t>l,^h) 


^ is-^^vA^a, 


(4.5b) 


as well as the linear map 






n 3 m;. ^ 


{mLC/'"^ox}. 


(4.5c) 



With the given bases for the spaces (4.3), the bilinear forms produce the matrix V^, 



and P+W^ + P~W^, while the linear form yields the vector /3o- 

Look around quadrature. What is missing to get a complete discrete set is the full 
discretization of the following bilinear forms 



and the linear form 



Tt 


X Sh3 


{l^l 


i'h) 


SI 


y-Th3 


ir. 


,Vh) 


SI 


xSh3 


(K 


i^h) 


SI 


^K 


1 — > 


{ 



{/^;„Ki/;,,}, 



(4.6a) 
(4.6b) 
(4.6c) 



{0;;,(V[/'-ox).n}. (4.6d) 

The elements of the space 5^*, can be decomposed as sums of elements of the spaces 



± 



5";^ := spanjxi : i e Zn} 



Therefore, the practical computation of all elements in (4.6) can be done if we are 
able to compute integrals of the form 



a+h/2 



{{t) ■ n{t)dt, 



a-h/2 



a+h/2 pb+h/2 



m{t, t) n(t) • n(T)dMT. 



(4.7) 



~h/2 Jb-h/2 



The approximation of integrals in one variable will be carried out with a three-point 



formula of order four using points outside the integration interval (see (3.2|) 



a+/i/2 



a-h/2 



(t){t)At 



h 
24 



{(l){a -h) + 22(t>{a) + 4>{a + h)). 



(4.8) 



Second collection of discrete elements. As already mentioned, the semidiscrete 



(4.7) 



elements (4.6) can be fully discretized once we approximate all integrals of the form 



For the one variable integrals we use (4.7) and for the double integrals we 

in each variable. Note that 



use the nine-point formula that arises from using (4. 



the normal vector appears always in the integration variable and that we have defined 
YVi := h n(tj), etc, which means that the value h will not appear in any of the resulting 
expressions. It is then easy to verify that this integration process transforms the 



bilinear forms (4.6a)~(4.6c|) into fully discrete bilinear forms associated to the matrices 

K;.,J;,andQ(M^,,- 

(4.6d) leads to the vector /3i 



Vj^ h)Q respectively. Finally, quadrature on the linear form 
in (|3^. 



5. Discussion on parameters. There are several choices related to parameters 
that we next proceed to discuss. The first parameter is the ±1/6 value that defines 
the staggered grids where the spaces 5'^ and T^ = span{(5j : i S Zat} are defined. 
These were first discovered in [3S] as the optimal choice of the parameter e such that 
the fully discrete method 



N 



log |x(^i - e/h)-yL{tj)\ Xj = g(ti -eh) i = 1, . 



,N 



provides a second order approximation of the parametrized Symm's equation 



log|x(T)-x(i)|A(t)di = 5(T). 



All other choices yield methods of order one, except e = which is not practicable and 
e — 1/2 which gives an unstable method. (Note that e + 1 leads to the same method 
as e.) With different techniques, these optimal choices were rediscovered in |S], where 
the method was shown to work for more complicated logarithmic kernels (such as the 
one for the Helmholtz equation), and where it was shown that e — ±1/6 were the 
only two values that led to second order methods. In fact, after some simplification, 
the leading term of the expansion in [Sj Proposition 16] is formally the quadrature 
error (the expansion holds in some Sobolev norm) 






N 



log#(i-T)u(T)dr- 



^/i^log_^(f- 



tj)u{tj) :^ hC {\og4 + \og^it/h)) u{t) + 0{h^) 



in terms of the periodic logarithmic function logjt(i) :— log(sin^(7ri)). (Note that this 
was also studied in ^ , where the logjt in the first order coefficient was not identified, 
although its graph was given.) This shows clearly that the best observation points for 
this quadrature error are those canceling the order one coefficient, namely, the points 
t = ihzL ^h, which are exactly the points that are used in our fully discrete methods. 
Only very recently T2], it was discovered (by the authors of the current paper), that 
the same structure could be used to find a Nystrom discretization of the hypersingular 
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operator written in integrodifferential form. In its turn, this led to the construction 
of two fuUy discrete Calderon Calculus of order two (one for each of e = ±1/6) in 

The values {jk, -rg) of the matrix Q come from the look-around quadrature formula 



(4.8). The need for using points around the integration interval in quadratures is re- 
lated to asymptotic behavior of the discretization errors: we want to have quadrature 
of sufhciently high order, but we do not want to introduce any more relative distances 
between points, since they would trigger first order asymptotic errors through the 



function C(e) := log 4 + log /^(e). We believe that the formula (4.8) might be new, but 
it has to be said that it has been derived in the same spirit as formulas in (TU] and 
[TT] . trying to keep fixed relative distances between integration points at the price of 
using points outside the integration interval. 



The following set of parameters is given by the definition of the fork (4.1) and 



the ziggurat (4.2), that is, they correspond to the matrices P^{a). The choice a — - 



6' 

was first discovered in [T^, applied just to the single layer operator V. The choice of 
parameters is motivated by the figure of the Dirac deltas fitting in a triangular shape 
(a basis function for the space of continuous piecewise linear functions). This is due 
to the origin of the method based on a variant of the qualocation methods of Ian 
Sloan |26j . In particular, the stability analysis for the corresponding matrix V/i (in 
form of an inf-sup condition [T31 Proposition 10]) is essentially outsourced to the work 
of Chandler and Sloan on qualocation methods [7] . A nice feature of the particular 
Dirac fork a — ^, following the shape of a hat function, is that its antiderivatives 
have the shape of the ziggurat, which mimics that shape of a B-spline of degree two, 
as corresponds to antiderivatives of hat functions. 

The interactions of Dirac deltas and characteristic functions can be expressed 
either with simple elements 

{xf,St) ■= 7, {x7+l,S^) = {xt-i,Si) := 1 - 7, {Sr,xf) := 0, otherwise, (5.1a) 

{^f,Xr)-^l, {St-i,Xi) = {Si+i,Xi) -=^-1, {Sf,Xi)-=0, otherwise, (5.1b) 

or with the composite actions of forks over simple characteristics and ziggurats over 
simple deltas (that is, with the elements of the mass matrix M): 

{S:,x^) = {xlS,) = l-2p, (5.2a) 

(54i,x.) = (x*±i,<5.)=p, (5.2b) 

{S*,Xj) = {x*i,Sj) = 0, otherwise. (5.2c) 



It is clear that, given the parameter a in (4.1|-(4.2), 7 determines p and vice versa 



What is less obvious, and we will try to explain next, is that a (and the choice of 
the quadrature rule), actually determines both sets of coefficients: p = jg^{7 — 6a) 
and 7 = ]^(1 + 12a). We start this argument with a simple computation. Let Tc be 
the translation operator T^X := A(- — c) and consider two collections of formal finite 
difference operators 

4^ := ^(T_j, + T5,) + ^{T_^.+T^^), Af := /3(r_, + T,) + (1 - 2/3)To. 



With this notation we can write 6* — A'^Si, x* — ^hXi and (4.8) becomes 



pa+h/2 

{X(a-h/2,a+h/2) , 0} = / (l>{t)dt « ^ ((/)(a - h) + 220(a) + <j){a + h)) 

J a-h/2 

= h{5a. ^TU] - /HAf '<5a, 0}. (5.3) 
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Similarly, the action of the composite difference operator A^A,/ = A^' A^ is 
given by the expression 

3L((21a+l)(r_i,,+ri,J + (22-21a)(T_|,+T|,J+a(r_|,,+rr,J + (f-«)(T_ii,+T^,j) 



(recall (3.4)). A simple computation then shows that 

Al/'%U-Kv = Oih*)^ ^ p=±{7-6a). (5.4) 

Let us now try to justify why ( |5.4[ ) is relevant. Imagine that we want to solve the 
trivial equation A = d^U™'^ with our class of methods. The non-conforming Petrov- 
Galerkin approximation of this equation is 

Ken, (x*,A^)=p(A,_i + A,+i) + (l-2p)A, = {x*,a„C/'"^} Vi. (5.5) 

The fully discrete method consists of separating x* into its ± parts and then using 
quadrature on each side. This leads to the following argument (see (5.3|): 

{X*,9„t/'"^} = {A^X.,5nf/'"^} = {%., A^a„C/'-} « /j{<5., A;^/'^A^5„C/'-}. 

The fully discrete realization of A = ^nC^'"*^ is then given by 

K^Y.^0^0^ Pi^^-l+K+l) + {l-'2p)K = h{S,,Al/^%'^^nU'^''}, Vz. (5.6) 



A dimensional look at (5.6) shows how the unknowns \j are trying to approximate 
hX(tj). The consistency error for equations (5.6) is then obtained when plugging in 
hX = hdnU"^'^ in the left hand side of the discrete equations and subtracting the 
right-hand side: h{{Al^X){U) - (A^/^'*A;^A)(i,)). This takes us back to ( [5^ . 

6. Building equations using the discrete calculus. We show here how 
to write integral equations for boundary value problems associated to the exterior 
Helmholtz equation: 

AU + k^U^O mn+, drU -I kU ^o{r-^/'^) at infimty. 

All formulations will be given directly at the discrete level. Here fi-j. is the exterior 
of a collection of smooth closed curves with non-intersecting interiors. 

Dirichlet problem. With a boundary condition ^U + 7[/'"'^ = 0, we can try four 
different formulations. In all cases, the trace of the incident wave is tested using (3.5 ). 
In the indirect formulations we have to give the integral equation and the potential 
representation. A single layer potential leads to an integral equation of the first kind 

YhV = (3o and Uh = Shi-)v, (6.1) 

while a double layer potential leads to an integral equation of the second kind 

^MiP + Kh'ip = f3a and Uh = Bhi-)ip. (6.2) 

In the direct formulations, we have a representation formula in terms of discrete 
Cauchy data: 

Uh^Sh{-)X-Dh{-)^. (6.3) 
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Here A can be found using one of two integral equations and ip will be derived by 
projecting data. We can use an integral equation of the first kind 

V,,A = -iMv3 + Kft¥', where Mtp = /3o, (6.4) 

or an equation of the second kind 

iMA + J/,A = -W;j¥', where My? = ,3o. (6.5) 

In both cases, A^ « VC/(mi) • n^. 

Neumann problem. Consider now a boundary condition d^U + d^U™'^ = 0, and 
test the incident wave as in (3.5) to produce a vector /3i. There are two possible 
indirect formulations: with the single layer potential 

- iMj7 + JftJ7 = /3i and f/,. = S,,( • )i7 (6.6) 

and with the double layer potential 

W;,^--/3i and C/^ = D;,(.)'0. (6.7) 



The direct formulations use the representation formula (6.31 and either the equations 
-\Mip + Kh^ = \h\ where MA = /3i, (6.8) 

or 

W/i(p = -iMA- J,,A, where MA = /3i. (6.9) 

In the direct representation ipi « ^U{m.i). 

Comhined potentials. If — fc^ is a Dirichlet eigenvalue of the Laplace operator in 
the interior domain n_ , then equations (6.1), (6.4), (6.6) and (6.8) are approximations 
of not uniquely solvable problems. Similarly, if — fc^ is a Neumann eigenvalue, all other 
four equations break down. Well posed equations for all frequencies can be found using 
a combined field integral representation: 

Uh = {'Dh{-)-ikSh{-))il, (6.10) 

leading to 

\Mr] + Khr]-ikYhri^lia (6.11) 

for the Dirichlet problem, and 

-^h'n + ^k\M■q-^kih'n^ 1^1 (6-12) 

for the Neumann problem. Direct formulations based on combined field equations can 
also be derived using the arguments of the Burton-Miller integral equation. 

7. Experiments in the frequency domain. Let Fi be parametrized by 

*^ (lI)'^) + 72((l+cos'(2^t))cos(2^t),(l+sin2(27ri))sin(27rt)) (\ " M , (7.1) 
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and let r2 be the ellipse parametrized by i i~> (4, 5) + (cos(27rt), 2 sin(27rt)). Discretiza- 
tion will be led by a single parameter N: we will take 2N points on Fi and N points 
on T2. We fix the wave number fc = 3 and consider a source point solution 



f7(x) 



i7^^^(fc|x-xo|) with 



xo 



vlO' 



10^ 



(7.2) 



Since the point Xq is in the interior of Fi, using U"^'^ = —U a.s incident wave, will 
give U as exact solution of the corresponding exterior problem. The boundaries of 
the scatterers are thus acting as transparent screens. We will measure errors 

E5,^*:= max|C/(z)-C/^(z)|, Obs = {(0,4), (4, 0), (-4, 2), (2, -4)}. (7.3) 

zGObs 

For direct methods involving the computation of A, we will compute 

E^ := N max|Aj - VC/(mj) • nj\. 
j 

The rescaling factor N is due to the fact that |nj| is proportional to h, instead of 
being of order one. For direct methods involving ip, we will compute 

E^ := max \(j)j — U{inj)\, where 4> = Qv'- 
j 



Note that the effective approximation of the trace in the discrete potential (3.1b) is 



not ip but (p — Q</3, which justifies our choice for the latter to compute norms of errors. 
It is clear that E™* measures the error of a smoothing postprocess and, as such, will 
benefit from weak superconvergence properties. On the other hand, the errors for 
the quantities on the boundary are measured in uniform norm. We will show that in 
all the experiments and for all the quantities, the errors are 0{N~^). Experimental 
orders of convergence are computed using errors on two consecutive meshes. 

First round of experiments. We first test all the formulations of Section [6] using 
the above geometry and exact solution. In all of them we test the simplest method 
(a = 1) and the method that generalizes the fork distribution in [T^ (a = 5/6), for 
which there is partial theoretical justification. Tables 7.1 to 7.10 show convergence 
of order three in all measurable errors. Note that the method for a = 5/6 is almost 
invariably slightly better than the method for a = 1. The errors are displayed in 
Tables |7.1| to |7.10[ corresponding to the ten integral equations given in Section [6j 



N 



a = 5/6 



e.c.r 



a = 1 



e.c.r 



10 


2.0504^(- 


-001) 




2.1788£:(- 


-001) 




20 


4.2900^(- 


-003) 


5.5788 


6.8665£:(- 


-003) 


4.9879 


40 


4.2678^(- 


-004) 


3.3294 


7.6927£:(- 


-004) 


3.1580 


80 


5.0466S(- 


-005) 


3.0801 


9.3497£;(- 


-005) 


3.0405 


160 


6.2217S(- 


-006) 


3.0199 


1.1603£;(- 


-005) 


3.0104 


320 


7.7503S(- 


-007) 


3.0050 


1.4477£;(- 


-006) 


3.0027 


640 


9.6795S(- 


-008) 


3.0012 


1.8087£;(- 


-007) 


3.0007 



Errors E^* for equation 



Table 7.1 
||6.1[l (indirect, single layer, Dirichlet). 



Tests on condition numbers. Equations associated to weakly singular and hyper- 
singular operators will have naturally growing condition numbers. In Figure |7.1| we 
show how cond(W/i) = 0{N), but cond(V/iW/i) = C(l), that is, the Calderon pre- 
conditioner works at the discrete level. We also show how integral equations of the 
second kind are well conditioned, by showing how cond(^M — K^) = 0(1). 
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N 


a = 5/6 


e.c.r 


a = ] 


L 


e.c.r 


10 


1.0885^(-001) 




1.1905£:(- 


-001) 




20 


2.1132^(-004) 


9.0086 


6.1488£:(- 


-004) 


7.5971 


40 


1.4713^(-005) 


3.8443 


4.0943£:(- 


-005) 


3.9086 


80 


1.5695^(-006) 


3.2288 


3.1627£;(- 


-006) 


3.6944 


160 


1.8971^(-007) 


3.0484 


2.8519S(- 


-007) 


3.4712 


320 


2.3782^(-008) 


2.9959 


2.9196i;(- 


-008) 


3.2881 


640 


2.9942S(-009) 


2.9896 


3.2775£;(- 


-009) 


3.1551 



Errors E™' for equation 



Table 7.2 
\6.2\ (indirect, double layer, Dirichlet). 



N 


a ===5/6 


e.c.r 


a = 


1 


e.c.r 


10 


1.1492^( 


-001) 




1.2300£:( 


-001) 




20 


9.5390^( 


-004) 


6.9125 


1.1919£:( 


-003) 


6.6892 


40 


1.1902^( 


-004) 


3.0026 


1.3916£:( 


-004) 


3.0985 


80 


1.4778^( 


-005) 


3.0097 


1.7282£;( 


-005) 


3.0095 


160 


1.8395^( 


-006) 


3.0060 


2.1610£;( 


-006) 


2.9995 


320 


2.2948S( 


-007) 


3.0029 


2.7039£;( 


-007) 


2.9986 


640 


2.8657^( 


-008) 


3.0014 


3,3822£;( 


-008) 


2.9990 


N 


a = 5/6 


e.c.r 


a — 


1 


e.c.r 


10 


4.5613S( 


-HOOO) 




4.6869£:( 


-HOOO) 




20 


2.3802S( 


-001) 


4.2603 


3.9297£;( 


-001) 


3.5761 


40 


1.9732^( 


-002) 


3.5925 


4.3200£;( 


-002) 


3.1853 


80 


2.3458^( 


-003) 


3.0724 


5.3704£:( 


-003) 


3.0079 


160 


2.8639^( 


-004) 


3.0340 


6.6578£:( 


-004) 


3.0119 


320 


3.5581^( 


-005) 


3.0088 


8.3179£:( 


-005) 


3.0007 


640 


4.4405^( 


-006) 


3.0023 


1.0395£:( 


-005) 


3.0003 



Table 7.3 
with exterior solution computed using | |6.3| l (direct, 
weakly singular integral equation, Dirichlet) . The upper table corresponds to E^' and the lower 
table corresponds to Ej^. 



Errors E™* 



and Ej^ for equation (|6.4|l 



Dependence with respect to a. It is unclear from the experiments whether there 
is a much better choice of the parameter a, that dictates the mixture of test functions 
in the method. Let us first show that a — 1/2 is not feasible. For a test equation 

as A'' increases. The domain is the curve 

It is clear from Figure 

However, 



(6.4) we compute the errors E 



and E5J?* 



Fi and the exact solution of the Helmholtz equation is (7.2) 



7.2 that E;^ is not converging, while E™* converges with the right order, 
inspection of the condition numbers show that they are of the order 10^'^. This makes 
the method highly unstable. Convergence of the potential solution can be explained 
by the fact that the potential postprocessing is a smoothing operator which, in some 
way, eliminates high frequency unstable components of the error and only observes 
approximation properties. In Figure |7.3[ we explore how the condition numbers of 
Vh blow up as a —> 1/2 and stay large (but considerably smaller) beyond this value. 

8. More complicated problems. 

8.1. Transmission problems. Consider now the domain fi interior to the curve 



(7.1) 



In addition to the exterior Helmholtz equation (2.3 1, we consider an interior 
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N 


a = 5/6 


e.c.r 


a — 


1 


e.c.r 


10 


3.4080^(-001) 




3.6686£:( 


-001) 




20 


l,9862^(-002) 


4.1008 


2.1362£:( 


-002) 


4.1021 


40 


1.2691^(-003) 


3.9682 


1.4680£:( 


-003) 


3.8631 


80 


8.3324^(-005) 


3.9289 


1.0955£;( 


-004) 


3.7441 


160 


6.1749^(-006) 


3.7542 


1.2007£;( 


-005) 


3.1896 


320 


5.7108^(-007) 


3.4347 


1.4394£;( 


-006) 


3.0603 


640 


6.7160-E;(-008) 


3.0880 


1.7631£;( 


-007) 


3.0293 


N 


a = 5/6 


e.c.r 


a — 


1 


e.c.r 


10 


2.5186^(+001) 




2.8815£;( 


+001) 




20 


1.5341^(+000) 


4.0371 


1.7131£;( 


+000) 


4.0722 


40 


1.2257^(-001) 


3.6457 


1.5322£:( 


-001) 


3.4830 


80 


9.2022^(-003) 


3.7355 


1.3286£:( 


-002) 


3.5276 


160 


7.9479^(-004) 


3.5355 


1.4155£:( 


-003) 


3.2305 


320 


7.9863^(-005) 


3.3150 


1.6693£:( 


-004) 


3.0840 


640 


9.3918^(-006) 


3.0880 


2.0384£:( 


-005) 


3.0338 



Errors E^' 



and EVf for equation 



Table 7.4 
l|6.5|l with exterior solution computed with 



6.3|l (direct, 



second kind integral equation, Dirichlet). The upper table corresponds to E^^* and the lower table 
corresponds to Ej^. 



N 



= 5/6 



e.c.r 



a = 1 



e.c.r 



10 


2.0581S(- 


-001) 




2.2349£;(- 


-001) 




20 


8.9154S(- 


-005) 


1.1173 


2.6185£;(- 


-004) 


9.7372 


40 


1.3500S(- 


-005) 


2.7233 


1.7638£;(- 


-005) 


3.8920 


80 


1.4812S(- 


-006) 


3.1881 


1.3681£;(- 


-006) 


3.6885 


160 


1.6462^(- 


-007) 


3.1696 


1.3916£:(- 


-007) 


3.2974 


320 


1.9224^(- 


-008) 


3.0981 


1.6450£:(- 


-008) 


3.0806 


640 


2.3189^(- 


-009) 


3.0514 


2.1447£:(- 


-009) 


2.9392 



Table 7.5 

61 (indirect, single layer, Neumann). 



equation with a different wave speed 

AV^ + (fc/c) V = in n. 
An incident wave C/'"'^ is given and two transmission conditions are imposed on F: 

j+u + Po^i~v, a+c/ + /3i = «a„F. 

In practical problems (/3o,/3i) := {'-fU™'^,dnU™'^). We clioose these transmission data 



so that the exact solution is the pair given by C/ in (7.2) and 

V{z):=eM<k/c)z-d), d:=(^,-^). 

We take k = 3, c ~ 2/3, and n — 3/2. The direct symmetric boundary integral 
formulation of Costabel and Stephan [5] is used. The unknowns are the Cauchy data 
for the interior problem, so that the integral representations are 



U = -S(fc)(A- - /3i) + D{k)iif- - /3o), 
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V^K-'S{I)X--B{^) 



v> 



N 


a ==5/6 


e.c.r 


a^ ] 


L 


e.c.r 


10 


1.4507^(-001) 




1.3385£:(- 


-001) 




20 


1.8995^(-002) 


2.9330 


1.9220i;(- 


-002) 


2.7999 


40 


9.3066^(-004) 


4.3512 


9.3855£;(- 


-004) 


4.3560 


80 


6.1122S(-005) 


3.9285 


6.1330£;(- 


-005) 


3.9358 


160 


4.3175S(-006) 


3.8234 


4.3356£;(- 


-006) 


3.8223 


320 


3.3804^(-007) 


3.6749 


4.2660£;(- 


-007) 


3.3453 


640 


3.0335S(-008) 


3.4782 


5.1771£;(- 


-008) 


3.0427 



Errors E^* for equation 



Table 7.6 
\6.7\ (indirect, double layer, Neumann). 



N 


a = 5/6 


e.c.r 


a = 


1 


e.c.r 


10 


1.8360^( 


-001) 




2.1658£:( 


-001) 




20 


3.2147^( 


-003) 


5.8357 


5.4219£:( 


-003) 


5.3199 


40 


3.2038^( 


-004) 


3.3268 


5.9516£:( 


-004) 


3.1874 


80 


3.7952^( 


-005) 


3.0775 


7.2500£:( 


-005) 


3.0372 


160 


4.6748^( 


-006) 


3.0212 


9.0125£:( 


-006) 


3.0080 


320 


5.8184S( 


-007) 


3.0062 


1.1255£:( 


-006) 


3.0014 


640 


7.2629S( 


-008) 


3.0020 


1.4067£;( 


-007) 


3.0001 


N 


a = 5/6 


e.c.r 


a = 


1 


e.c.r 


10 


3.3579S( 


-001) 




3.6830E{ 


-001) 




20 


9.7882S( 


-003) 


5.1004 


1.6541£;( 


-002) 


4.4768 


40 


9.8787^( 


-004) 


3.3087 


1.9671£;( 


-003) 


3.0719 


80 


1.1104^( 


-004) 


3.1533 


2.4081£;( 


-004) 


3.0301 


160 


1.3330^( 


-005) 


3.0583 


3.0099£;( 


-005) 


3.0001 


320 


1.6404^( 


-006) 


3.0226 


3.7716£:( 


-006) 


2.9946 


640 


2.0374^( 


-007) 


3.0092 


4.7276£:( 


-007) 


2.9960 



Errors E^^' 
integral equation, Neumann). The upper 

to Etr. 



Table 7.7 

Si, with potential representation ( |6.3| (direct, second kind 
table corresponds to E™* and the lower table corresponds 



The corresponding system of integral equations is 



W{k) + KW{^] 
-K(A;)-K(f) 



J(fc) 
V{k) + 



-J(f) 






W(fc) 



^I + J(A:) 



il-K(fc) V(fc) 



Pi 



We discretize each of the elements in the system of integral equations and in the 
integral representations using the rules of the discrete Calderon Calculus. Taking TV 



discretization points on the boundary, we compute the exterior error (7.3) and errors 
on the boundary 



E^ := N max \Xj - k VF(mj) • n^ 



^N ■ — 



max I (j)j 
j 



Vim, 



with (p = Qip. 



The corresponding errors are plotted in Figure 8.1 



8.2. CQ discretization in the time domain. In this final example we show 
how to combine the fully discrete Calderon Calculus with a Convolution Quadrature 
routine to produce time-domain discretization of scattering of waves by obstacles. We 
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N 


a ==5/6 


e.c.r 


a — 


1 


e.c.r 


10 


1.7474^( 


-001) 




1.5968£:( 


-001) 




20 


7.0648^( 


-003) 


4.6284 


8.6420S( 


-003) 


4.2077 


40 


4.5988^( 


-004) 


3.9413 


6.3940£;( 


-004) 


3.7566 


80 


4.2497^( 


-005) 


3.4358 


6.5G02£;( 


-005) 


3.2982 


160 


4.4521^( 


-006) 


3.2548 


7.2760£;( 


-006) 


3.1593 


320 


5.0449^( 


-007) 


3.1416 


8.5836£;( 


-007) 


3.0835 


640 


5.9863£'( 


-008) 


3.0751 


1.0416£;( 


-007) 


3.0428 


N 


a = 5/6 


e.c.r 


a — 


1 


e.c.r 


10 


1.0240^( 


+000) 




9.9G18£;( 


-001) 




20 


1.0248^( 


-001) 


3.3207 


1.1G62£;( 


-001) 


3.1621 


40 


5.8839^( 


-003) 


4.1225 


7.3167£:( 


-003) 


3.9182 


80 


4.3947^( 


-004) 


3.7429 


6.3631£:( 


-004) 


3.5234 


160 


3.7752^( 


-005) 


3.5411 


6.3618£:( 


-005) 


3.3222 


320 


3.7134^( 


-006) 


3.5457 


7.0366£:( 


-006) 


2.1765 


640 


4.0517^( 


-007) 


3.1962 


8.2524S( 


-007) 


3.0920 



Errors E™* and E^ for equation | |6.9| l 



equation, Neumann). 



Table 7.i 
with potential representation 



6.3| (direct, hypersinqular 



The upper table corresponds to E^' and the lower table corresponds to E|^. 



N 


a = 5/6 


e.c.r 


a = ] 


L 


e.c.r 


10 


1.2545^(-001) 




1.3462£;(- 


-001) 




20 


1.4698^(-003) 


6.4153 


2.1018£:(- 


-003) 


6.0011 


40 


2.7686^(-004) 


2.4084 


4.6628£:(- 


-004) 


2.1723 


80 


4.2480^(-005) 


2.7043 


7.6357£:(- 


-005) 


2.6104 


160 


5.7877^(-006) 


2.8757 


1.0648£:(- 


-005) 


2.8422 


320 


7.5027^(-007) 


2.9475 


1.3925£:(- 


-006) 


2.9348 


640 


9.5326^(-008) 


2.9765 


1.7758£;(- 


-007) 


2.9711 



Table 7.9 
for equation l |6.11| with potential representatii 
potential, Dirichlet). 



Errors EJ^' 



(|6.10| (indirect, combined field 



first explain some general ideas of the CQ method. More details, specifically applied 
to scattering problems, are given in [2li21|. while the original ideas of multistep-based 
CQ for hyperbolic problems appear in [22]. 

Generalities about CQ. We start with a causal approximation of the derivative: 
if K > 0, then the operator 



9. 



^{lu~ui--K) + lu{--2K)) 



(8.1) 



is the backward differentiation operator associated to the BDF2 method. The associ- 
ated transfer function (the Laplace transform of the operator) is 



H 



le" 



-2k,s 



(8.2) 



Let now Afi{s) be any of the elements of the discrete Calculus (one of the potentials 
or one of the operators), with k = —is, s € C and Res > 0. This is the same as 
saying that we are taking the operators associated to the Laplace resolvent equation 



AU- 



'-U 



in M^\r (radiation conditions are reduced to imposing U G _ff^(M^\r), 
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N 



a = 5/6 



a^ 1 



10 1.1559^(-001) 1.2167£:(-001) 

20 2.7343^(-003) 5.4017 3.4492£:(~003) 5.1406 

40 9.5433^(-005) 4.8405 1.3958£:(-004) 4.6271 

80 5.6897^(-006) 4.0681 9.9258£:(-006) 3.8138 

160 4.0825^(-007) 3.8008 6.9688S(-007) 3.8322 

320 3.2007^(-008) 3.6730 5.0311S(-008) 3.7920 

640 2.9313^(-009) 3.4488 4.0423£;(-009) 3.6376 
Table 7.10 

Errors E™' for equation |6.12| with potential representation ( |6.10| (indirect, combined field 
potential, Neumann). 




condllion(1/2M-K) 



450 500 



150 2aa 250 300 



Figure 7.1. Condition numbers for the matrices Wh, V^Wh and ^M — K;,. The results are 
given for the choice a = 5/6. Results for a = I are almost identical. 

which in practice imposes exponential decay at infinity). After some manipulation in 
the complex plane, we can write 



^hisf,) = ^ A^^h[m]e 



m=0 



The Convolution Quadrature method is the practical computation of convolutions of 
the form 



A,,(9k)i/) = ^ Ak,/iH'/'(- - "ik) 



(8.3) 



(compare with (8.1) and (8.2)). The forward convolution form consists of sampling a 
causal function t/; : M — >■ C^^ , denoting ^p[n] := xp{Kn), and then computing 



Ah{d^)ip[n] := ^ A^ji[m]xp[n - m]. 



(8.4) 



T7l — 



(Note that we use the same notation, but now ip is discrete in time, i.e., it is a sequence 
of vectors.) The same idea can be used to solve convolution equations (in the same 
way that (lO) is the seed of the BDF2 method) 



^^ AKji[m]xp[n - m] ^ $,[n] n = 0, 1,. 
19 



(8.5) 



m=0 





Figure 7.2. The figure on the left shows history of convergence for the choice a = 1/2 using a 
direct single layer potential based method. The method is clearly not converging for the unknown on 
the boundary, but convergence is restored in the smoothing postprocessing of the potential. The figure 
on the right shows a history of convergence w.r.t. a for fixed N. The peak at a = 1/2 corresponds 
to the unstable choice of this parameter. 





Figure 7.3. Condition number of the matrix V^ as a function of the parameter a. The choice 
a = 1/2 equalizes the height of the four Dirac deltas in Figure \4-.1\ making the method unstable. 
Past this threshold, condition numbers are unreasonably high. 



where ^ : M — > C^ is a given causal function sampled at the points Kn, or $\n\ are 
the entries of a sequence of vectors ^. Note that in (8.4 1 and (8.5) data (■0 and 
^ respectively) are sampled in the time domain, while the action of the operator is 
taken using the transfer function. Practical ways of computing these convolutions are 
explained in [2]. They involve a clever use of FFT, contour integrals, and multiple 
evaluations of the transfer function Kh{s). In the case of the convolution equation 



.5), repeated inversion of A^ /i[0] 



Aft(|^) is also required. 



A scattering problem. In this first example, we use the time domain version of 
(6.3) and (6.9). The normal derivative of an incident plane wave C/'"'^(i,x) is sampled 



at the observation points at all times 

(3f [n] := -(VC/'"'=(n k, m±) • rii, . . . , V?7'"'=(n k, m± ) • Ha,)^, n > 0. 

We assume that the discrete function /3i[n] :— P+/3j'"[n] + P^/3f [n] is causal: this is 
true in the reasonable physical situation when the incident wave has not reached any 
of the obstacles at time zero. We then solve equations looking for causal sequences 
if = {ip[n\) and A = (A[n]) satisfying 

M\[n] - /3i[n], WhidkMn] = -iMA[n] - J,,(a«)A[n], Vn > 0. (8.6) 

20 





Figure 8.1. Errors for the transmission problem. On the left, errors for the choice a = 1. On 
the right, for a = 5/6. 



The potentials are then computed at every time step using the CQ method once again, 
resulting in sequences 

U[n] = Sh{d^)X[n] - Dh{d^)cp[n]. (8.7) 

Note that this is a fully discrete method for the scattering of a sound-hard obstacle 



by a transient incident wave. Note also that the sequence of functions (8.7) are a 
classical solution of the BDF2-discretized wave equation [22] : 



dlU[n] - AU[n] =0 in M^ \ r 



Vn. 



To test the method, we change some signs so that we end up solving an interior 
boundary value problem, namely, we solve W^(9k)¥' = ^ MA — 3h{d^)\, instead of 
the second equation in (8.6). The potential solution (8.7) is then an approximation 
of-C/'"=(nK, •) inrj_. 

For the experiments we take the boundary of the domain parametrized with 

^(4 (1 + cos2(27rt)) cos(27ri), 5(1 + sin2(27rt)) sin(27rf)) (J T^M , 

the incident wave given by 

U"'%t,z):=p{t-R+z-d), R=1.2, d:=(-^,-^), p(t) := sin3(3i)xt>o, 

N discretization points on the curve, and M time steps of length T/AI, where T = 5. 
Finally we compute errors 



E% ^j := max |0,[M] + f/-^(r, m,)|, 



Zo = (0.2,0.2), 

cl)[M] = QifiM]. 



The values of N and M are chosen so that 0{N-^) = 0{M-'^): for j = 10, . . . , 19, 
we define 



iV, = L20(1.2)^"J, 



M, 



[iV3/2 20-i/2j^ Nf^20Mf 



N,+,/N, « 1.2. 



The results are reported in Table |8.1[ Experimental convergence rates are shown to 
confirm that the errors in 0{N^^). 
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N 


M 


TT^CXt 


e.c.r 




c.c.r 


123 


305 


7.1971£;( 


-002) 




1.3079£;( 


-001) 




148 


402 


4.2559£;( 


-002) 


2.8816 


7.6194£;( 


-002) 


2.9634 


178 


531 


2.4632£;( 


-002) 


2.9994 


4.3811£;( 


-002) 


3.0353 


213 


695 


1.4327£;( 


-002) 


2.9723 


2.5594£;( 


-002) 


2.9482 


256 


915 


8.2648£;( 


-003) 


3.0173 


1.4754£;( 


-002) 


3.0211 


308 


1208 


4.7404£;( 


-003) 


3.0489 


8.4560£;( 


-003) 


3.0532 


369 


1584 


2.7533£;( 


-003) 


2.9801 


4.9135£;( 


-003) 


2.9776 


443 


2084 


1.5894£;( 


-003) 


3.0135 


2.8368£;( 


-003) 


3.0129 


532 


2743 


9.1716£;( 


-004) 


3.0157 


1.6368£;( 


-003) 


3.0162 


638 


3603 


5.3072£;( 


-005) 


3.0005 


9.4818£;( 


-004) 


2.9946 



Errors E^^j^j and E^ j^j 



Table 8.1 
for an interior problem in the time dor 



A final experiment. To illustrate the capabilities of the time-domain discretiza- 
tion, we choose a kite-shaped sound-hard obstacle, hit by a short plane incident wave, 
and we plot several snapshots of the total wave field (incident plus computed wave). 
Results are shown in Figure |8?2) 
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t ^ 7.6560 





t = 8.3760 



t = 9.0960 





Figure 8.2. Six snapshots of the scattering of a plane wave by a kite-shaped sound-hard obstacle. 
The profile of the wave can be oberserved in the first two images, as it travels to the right and to the 
top. Discretization has been carried out with the order three Calderon Calculus and a BDF2-based 
Convolution Quadrature routine. 
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